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ABSTRACT 

Aims. Periodic variability in young stellar objects (YSOs) offers indirect evidence for an active dynamical mechanism. Starspots, 
accretion, stellar companions, and disk veiling can contribute to the photometric variability of YSOs. 

Methods. As part of an ongoing study of the p Oph star forming region, we report the discovery of 92.6 day periodic variations for 
the Class I YSO YLW 16A, observed over a period of three years. A SED model was fit to available photometric data for the object. 
Results. We propose a triple-system with an inner binary with a period of 93 days eclipsed by a warped circum-binary disk. The 
nature of the secondary is unconstrained and could be stellar or sub-stellar We report the discovery of a tertiary companion at a 
projected separation of ~40 AU that could account for the circum-binary disk warp. This light curve and model is similar to the model 
we proposed for WL 4 in previous work. Understanding these systems may lead to insights about the nature of stellar evolution and 
planetary formation, and provide valuable benchmarks for future theoretical modeling and near- and mid-infrared synoptic surveys of 
YSOs. 
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I 1. Introduction 

^ Star formation involves the gravitational collapse of a massive 

^ cloud core. Between this initial collapse, and the final con- 

I traction onto the main sequence, the protostar is classified as 

a young stellar object (YSO). These YSOs have ages of a few 

I— i million years (~1-10 Myr), and are characterized by high lev- 

^ els of accretion, ejection, and magnetic activity, as well as pho- 
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tometric variability (Joy 1945). The evolution of YSOs falls 
broadly into four stages. Class objects consist of a collaps- 
ing cloud core. Class I objects are protostars embedded inside 
a spherically-symmetric infall envelope. Class II objects, also 
known as "classical T Tauri stars," contain a stable primordial 
disk. The dispersion of the stable disk reveals a Class III object, 
or diskless "weak-lined T Tauri star" (e.g. Adams et al. 1987). 
. . The spectral energy distribution (SED) of YSOs differ from nor- 
^ mal stars by exhibiting an infrared excess as the circumsteilar 
l^ material reprocesses the central radiation. The amount of excess 
r> is strongly correlated with the evolutionary stage of the YSO. 
i-rt Due to the strength of the infrared emission, it is natural to study 
YSOs at these wavelengths. 

It is widely accepted that planetary systems form out of pri- 
mordial protostellar disks, and because such disks are an essen- 
tial structure in the evolution of Class II YSOs, the study of 
YSOs can lend valuable insights into the processes by which 
planets form (Lin & Papaloizou et al. 1980; Ida & Lin 2010, and 
references therein). At optical wavelengths, some YSOs are ob- 
served to exhibit periodic photometric variability (e.g., Rebull 
2001; Covey et al. 2006). The observed variability is generally 
attributed to the rotational modulation of large cold spots, hot 
spots, accretion and disk veiling. Photometric variability driven 
by rotational modulation of the proto-star are less pronounced 
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at infrared wavelengths, thus improving sensitivity to variabil- 
ity driven by disk-related processes and the subject of many re- 
cent NIR and Spitzer Space Telescope studies such as YSOVAR 
(Morales-Calderon et al. 2011; Flaherty et al. 2012; Flaherty & 
Muzerolle 2010; Faesi et al. 2012). 

p Ophiuchus (p Oph) is a nearby (~135 pc) star-forming 
region containing a few hundred such YSOs from a few Solar 
masses down to the free-floating planet mass regime (Mamajek 
2008; Marsh et al. 2010). In this paper we investigate the YSO 
in p Oph: YLW 16A. YLW 16A ( = IRAS 16244-2432, 2MASS 
J16272802-2439335, ISO-Oph 143, IRS 44) is classified as a 
Class I protostar (e.g. Luhman & Rieke 1999; Barsony et al. 
2005) which has been a notable subject of a previous study at X- 
ray wavelengths (Grosso 2001; Imanishi et al. 2001). Imanishi 
et al. (2001) detected an unusual bright X-ray flare, with a peak 
luminosity of 1.3x10"" ergs s"'. A 6.4 keV emission line was 
identified, which was attributed to fluorescence of cold neutral 
iron in the circumsteilar gas. An extended ~3400 AU (Beckford 
et al. 2008) nebulosity has been observed around YLW 16A in 
the infrared {H and K^ bands; Simon et al. 1987; Lucas & Roche 
1998) and at thermal radio wavelengths (Leousetal. 1991;Girart 
et al. 2004). High-resolution HST NICMOS imagery, obtained 
June 1998, reveals two nonpoint sources separated by 0.5", with 
flux ratios of 1.5 at 1.1 pva and 1.1 at 1.6 /im (Allen et al. 2002). 
Beckford et al. (2008) interpret the second source as being due 
to a reflection from a dusty jet inside a bipolar cavity. However, 
their conclusions do not rule out the possibility of a binary com- 
panion. Herczeg et al. (2011) also find a resolved binary with 
CRIRES/VLT, finding the west component to be 0.69±0.12 mag 
fainter at M-band, with CO and extended H2 emission, but no 
CO emission from the east component. Simon et al. (1987) re- 
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ported a variability of ~ 1 .0 mag in the K, band, over a timescale 
of approximately six months. Finally, Evans et al. (2009) derive 
an extinction of Ay = 9.8 mag for YLW 16A from the Spitzer 
c2d survey, which would correspond to a Aj~2. 

Doppmann et al. (2005) presented a high S/N near-infrared 
spectrum of YLW 16A, which is essentially featureless. There 
is, however a sharp H2 emission line at 2.1218 //m, characteris- 
tic of still accreting YSOs, and weak absorption features includ- 
ing the CO 2-0 bandhead at 2.3 jum. In general. Class I spectra 
tend to exhibit weaker features due to veiling and continuum 
emission characteristic of a protostar surrounded by a thick en- 
velope (Greene & Lada 2000). Covey et al. (2006) measure a 
local standard rest velocity of 4.41 km/s, and do not identify evi- 
dence for a double-lined spectroscopic binary companion, down 
to the ~l-2 km/s level, from the structure of the radial velocity 
cross-correlation function in a single epoch spectrum (Covey et 
al. 2013). 

Another YSO located in p Oph is WL 4, a Class II ob- 
ject, whose periodic photometric variability was discovered by 
Plavchan et al. (2008a) to be 130.87±0.40 days. The authors 
suggest a triple-YSO model, consisting of an inner binary and a 
third companion further out. A circum-binary disk, tilted with 
respect to the binary's orbital plane by the gravitational influ- 
ence of the third companion, eclipses each member of the inner 
binary in turn. Periodic eclipsing of a binary by a circum-binary 
disk is also the preferred model of the well-studied system KH- 
15D (Herbst et al. 2010, and references therein), as well as the 
recently discovered object CHS 7797 in the Orion star-forming 
region (Rodriguez-Ledesma et al. 2013, 2012). 

In this paper, we present the discovery of periodic near-IR 
photometric variability for YLW 16A. Our analysis parallels 
much of the analysis in Plavchan et al. (2008a). We present 
our observations in §2, and results in §3. In §4, we discuss the 
implications of our observations and proposed model. This dis- 
covery demonstrates that systems like YLW 16 A, WL 4, KH- 
15D and CHS 7797 may be common in multiple star- forming 
regions, constituting a new class of disk eclipsing YSOs. These 
systems will be valuable to study the evolution of circumstellar 
disks around YSOs, and potentially the formation sites for cir- 
cumbinary planets (Doyle et al. 201 1). 



2. Observations 

The photometry for the J, H, and K, bands were obtained from 
the Two Micron All-Sky Survey (2MASS) Calibration Point 
Source Working Database (Cal-PSWDB) (Skrutskie et al. 2006). 
Between 1997 and 2001, 2MASS imaged the entire sky in three 
near-infrared bands, J, H, and Kg. Hourly observations of 35 dif- 
ferent calibration fields were used to calibrate the 2MASS pho- 
tometry. One such field is located in p Oph, 8.5' wide in R. A. by 
60' long in deck, and centered at (R.A., decl.) = (246.80780°, 
-24.68901°). A total of 1582 independent observations were 
made of this field, which contains YLW 16 A, as discussed in 
further detail in Plavchan et al. (2008a); Parks et al. (2013). 

Two NACO images of YLW 16 A, taken from the European 
Southern Observatory (ESO) archive for the instrument, were 
obtained for the Ks and L bands (Lenzen et al. 2003; Rousset et 
al. 2003), as seen in Figure 1. The L band image was obtained 
on 9 April 2005, while the Ks band was obtained on 30 April 
2005. The NACO images allow us to measure a flux ratio be- 
tween YLW 16AA and YLW 16AB. Using aperture photometry, 
we derive flux ratios from the NACO images of 0.22 and 0.98 
at Ks and L respectively. We calibrate to the total flux from the 



2MASS Ks magnitude in the faint state (for the Ks NACO im- 
age) of 51.5 mJy and IRAC 3.6 microns in the faint state (for L) 
of 695.8 Jy. 

Photometry for YLW 16A was also obtained from the IRAC 
(3.6, 4.5, 5.8, and 8 //m) and MIPS (70 //m) instruments on the 
Spitzer Space Telescope, as part of the Cores to Disks (c2d) 
Spitzer Space Telescope Legacy program (Evans et al. 2003; 
Padgett et al. 2008). YSOVAR has also obtained photometric 
time-series of YLW 16A at 3.6 and 4.5 jum during the Spitzer 
Space Telescope warm mission, which will be part of a separate 
future publication (Morales-Calderon et al. 201 1). 

Finally, photometry was obtained from the literature at 10.8 
/im (Barsony et al. 2005), 850 /im (Jorgensen et al. 2008), and 
1.2 mm (Stanke et al. 2006). Excluding the NACO photometry, 
all other photometry is a blend of the YLW 16A system. The av- 
erage 2MASS photometry, as well as the additional photometry 
values are summarized in Table 1 . 



3. Analysis and Results 

3. 1 . Periodic Variability 

Periodic variability is readily apparent from a visual inspec- 
tion of the YLW 16A light curve (Figure 2). Using the Lomb- 
Scargle periodogram (Scargle 1982), the Box Least Squares pe- 
riodogram (Kovacs et al. 2002), as well as the period-searching 
algorithm of Plavchan et al. (2008b), a period of 92.62±0.84 
days is identified (also see Parks et al. (2013)). Access to all 
three algorithms are available online in interactive form at the 
NASA Exoplanet Archive, and include rudimentary estimates of 
the false-alarm probabilities (p-values). The period and period 
1-cr error are assessed from the Plavchan periodogram peak and 
Half- Width Half-Maximum respectively. All three algorithms 
detect the signal at 92.6 days, and period aliases of one-half and 
two times this period. 

The bin-less phase dispersion minimization Plavchan peri- 
odogram is adept at detecting arbitrarily-shaped periodic signals, 
compared to sinusoids for Lomb-Scargle and box-like transits 
for BLS, and thus detects the 93-day period with higher statisti- 
cal significance (5.4-cr). All three algorithms detect aliased pe- 
riods at integer fraction multiples of one day (e.g. periods of 
1/4,1/3,1/2,2,3, & 4 days,etc.) with statistical significance that 
in the case of Lomb-Scargle exceeds the 93-day signal. This is 
due to the long-term (e.g. ~400-500 day) variations seen in the 
light curve aliased with the ~ 1 observation per day cadence. Vi- 
sual inspection of the phased time-series confirms that these are 
aliased "false" periods. All three algorithms detect these "red- 
noise" long-term trends that are likely astrophysical in origin, 
but we do not quantify this time-scale further (Parks et al. 2013). 

Figure 3 shows two repeated cycles of the phased light and 
color curves of YLW 16A. For JD of 2,450,000.0 (note, not 
MJD=0), the corresponding phase is in Figure 3. The aver- 
age variation between bright and faint states is ~0.5 mag in the 
Ks band, with a maximum range of aKs - 0.95 mag. There is 
also periodic variability in the (H - K,) curve. The shape of the 
color variations diff'er from the shape of the photometric varia- 
tions in the Ks band, with an approximately sinusoidal phased 
curve with an amplitude of ~0. 15 mag and a maximum range of 
a(H - Ks) - 0.34 mag. The mean Ks magnitude is 10.22 and 
the mean (J - Kg) color is 7.07. The low S/N / band data is 
due to the high extinction towards the system. There are two 
data points that appear to be particular errant from the phased 
time-series in Figure 3 at phases of ~0.7 and 0.8, which may in- 
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dicate additional sporadic variability such as a flare that was not 
adequately sampled by the cadence of our observations. 

The folded light and color curves show that YLW 16A has 
both dim and bright states similar to WL 4 as discussed in 
Plavchan et al. (2008a). However, the bright state of YLW 16A 
is about half as long in phase duration as that of WL 4, such that 
it appears very much like an "upside-down eclipse". Addition- 
ally, the source is redder in // - K^ when the bright state occurs, 
which runs counter-intuitive to models for extinction variability 
(Parks etal. 2013). 

3.2. SED Modeling 

A model SED fit is generated, as done in Plavchan et al. (2008a). 
Our model SED can include up to three stellar components and 
five dust components, each with an independent extinction mag- 
nitude. We also allow for a variable broken power law for the 
system extinction as a function of wavelength (e.g. A^ oc ^«i'«2 
for A <,> /lo), and we dynamically color-correct 24 and 70 jum 
Spitzer photometry. Reddened Phoenix NextGen (Hauschildt 
et al. 1999) spectra contribute the stellar component(s) of the 
SED model, while the dust emission is modeled as a blackbody. 
By varying the temperatures of the stars and dust, as well as 
their effective radii, the composite curve was obtained using a 
X^ minimization fit to the available photometric data of both the 
bright and faint states. Given the multi-parameter fit, we have 
improved upon the analysis of Plavchan et al. (2008a) to include 
an AMOEBA simplex code parameter optimization (Nelder & 
Mead 1965). Figure 4 shows the best model SED fits to each 
of the bright and faint state SEDs, and the best fit parameters are 
given in Table 2. Simplex codes are generally sensitive to the 
initial guesses. Thus, our fit likely represents only one of multi- 
ple degenerate models that can adequately describe the SED in 
both the bright and faint states. However, we have chosen to fix 
some of the parameters to enforce ensure that these parameters 
remain consistent between the bright and faint states, and given 
the limited degrees of freedom. 



4. Discussion 

From the NACO images in Figure 1, we confirm the presence of 
at least two YSOs in the YLW 16A system. The L band NACO 
image shows two sources of approximately equal brightness, 
whose projected separation is around 0.3", corresponding to a 
projected separation of ~40 AU at a distance of 135 pc. The K^ 
band NACO image shows a more complex nebulosity surround- 
ing the visual binary, indicating at least one component of the bi- 
nary is deeply embedded in a protostellar envelope. Both NACO 
images were obtained in the faint state (L band phase =0.4656, 
Ks band phase =0.6920). These images do not constrain which 
of the visual components is the source of the observed photomet- 
ric variability. However, if the stellar components are of equal 
spectral type and radii, the fainter K-band component is a likely 
candidate given that the system was in a faint state at the time. 
Time-resolved adaptive optics monitoring is necessary to con- 
firm which component contributes to the photometric variability. 
Our SED model fit confirms the Class I nature of YLW 16 A. 
A fit with two dust blackbody components and no stellar com- 
ponents yields an unphysical extinction power law, and thus the 
SED model necessitates at least one stellar component. The large 
extinction of the stellar components and potential differences in 
infrared excess do not motivate/yield any useful constraints on 
multiple component stellar radii nor temperatures, even though 



the NACO images indicate at least two components. For exam- 
ple, a higher stellar temperature is partially degenerate with a 
larger stellar extinction. Rather, we derive a single composite 
set of stellar parameters from the fit of T* ~3500 K and effective 
stellar radius of ~4.66 Rq. If the system consists of three stars 
with equal radii and temperatures in the bright state, the corre- 
sponding stellar radii of each of the three components would be 
~2.7 Rq. If the system consists of two stars with equal radii in 
the faint state - e.g., a third star is completely extincted/eclipsed 
by a primordial disk - the corresponding stellar radii of each of 
the two components would be ~3.0 Rq. These radii are plausible 
for Class I YSOs. Our value of Ay ~ 10 yields a visual extinc- 
tion much larger than that derived in Evans et al. (2009), despite 
providing an initial guess for the stellar extinction in our model 
SED fit of Ay = 2. 

The 92.6 day periodic photometric variability cannot be as- 
sociated with the Keplerian orbit of the projected ~40 AU visual 
companion. Invoking the discussion in Plavchan et al. (2008a), 
the periodic photometric variability of YLW 16A is not readily 
attributable to starspots, chaotic disk extinction from accretion, 
nor other stellar activity induced variations that tend to operate 
on time-scales of less than a week. Instead, we postulate that 
the long-term periodic photometric variability of YLW 16A in- 
dicates the presence of a tertiary companion of unknown mass 
within a few AU of one of the visual components and with an 
orbital period of 92.6 days. If the tertiary companion is approx- 
imately the same mass/brightness of the other two companions, 
and it is periodically eclipsed by a circum-binary disk, this sce- 
nario could explain both the periodic variability as well as the 
~l/3 reduction in flux between bright and faint states. The lack 
of a detection of binarity in Covey et al. (2006); Doppmann et 
al. (2005) indicates that this tertiary companion may instead be 
at a much smaller observed luminosity. The SED model indi- 
cates that the luminosity of the hot dust component must also 
change substantially to explain the observed variability at IRAC 
mid-IR wavelengths. In other words, at K^ and IRAC bands, we 
are seeing possible periodic eclipses (shadowing) of some of the 
hot dust material in the system associated with this tertiary com- 
panion, rather than the proto-star photosphere itself. This in turn 
implies a strong star-disk dynamical interaction. 

The literature photometry for YLW 16A provides an indica- 
tion of the stability of the system, supporting a Keplerian ori- 
gin to produce the observed periodic variations. K^ band pho- 
tometry from Simon et al. (1987) and Wilking et al. (1989), as 
given in Table 1, indicates that the bright state magnitude of ~9.8 
may have been consistent for over twenty years. However, the 
faint state value of ~8.8 mag does not agree with our observa- 
tions. This could be due to long-term evolution in primordial 
disk structure. Follow-up observations with the Spitzer Space 
Telescope and the YSOVAR program will further constrain the 
long-term stability of the observed periodic variability (Morales- 
Calderon etal. 2011). 

Finally, the discovery of a system similar to WL 4 indicates 
that such systems may be common, and more may be uncovered 
with long-term photometric NIR and mid-IR intense photomet- 
ric monitoring. Both WL 4 and YLW 16A possess visual com- 
panions at wide separations. This implies that a wide companion 
offers a plausible, if not required, mechanism to explain how a 
circum-binary disk could be warped with respect to the orbit of a 
hypothesized inner binary, to produce the observed periodic disk 
extinction. These two systems join KH-15D in NGC 2264 and 
CHS 7797 in the Orion star- forming region (Herbst et al. 2010; 
Rodriguez-Ledesma et al. 2013, 2012). 
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5. Conclusion and Future Work 
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Fig. 1. NACO imagery of YLW 16A in the K, (Upper) and L (Lower) 
bands. The separation of the two sources in the L band is approximately 
0.3" (~40 AU projected separation). Both these images were obtained 
in the faint state (L band phase =0.4656, K band phase =0.6920). Color 
bars shown correspond to non-normalized counts. North is up and East 
is to the left. 
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Fig. 2. Top three panels: the J, H, & A'j-band light curves for 
YLW 16A, generated using data from the 2MASS Cal-PSWDB. "Scan 
groups" of six measurements taken in 10 minutes of elapsed real time 
are co-added, as in Plavchan et al. (2008b). Propagated l-cr error bars 
are shown in teal. Bottom tliree panels: Plavchan, Lomb-Scargle and 
Box Least Squares periodograms of the A'j-band light curve (Scargle 
1982; Kovacs et al. 2002; Plavchan et al. 2008b) 
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Fig. 3. From top to bottom: the /, H, K„ H-K„ & J-H phased light 
and color for YLW 16A, generated using data from the 2MASS Cal- 
PSWDB as in Figure 2, folded to a period of 92.6 days and plotted as a 
function of phase. A second phase of the same data is repeated. 
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Fig. 4. The model SED fit to observed YLW 16A photometry. Blue 
circles correspond to JHK and Spitzer IRAC photometry during the 
bright state; red circles correspond to the faint state. The summed model 
SED fit (black curve: bright state; maroon dashed curve: faint state) has 
contributions from a composite star (lower synthetic SEDs peaking at 
~0. 1 Jy; red: bright state; blue: faint state), hot dust (upper left black- 
body curves peaking at a few Jy, red: bright state; blue: faint state), and 
cold dust (green blackbody curve). Ground-based historical 10.8 ;um 
photometry, Spitzer MIPS 70 micron photometry and sub-mm photom- 
etry is shown as black circles (YLW 16A is saturated at 24 pm in the 
Spitzer c2d survey data). 
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Table 2. SED Model Parameters 

Parameter Value 

Fixed 



Distance 

short-zl extinction power law a 

\ong-A extinction power law a 


135 pc 

-2.0" 

-1.0" 


Varying in faint state fit, fixed in bright state fit to faint fit 


A Extinction transition 
cold dust T, L 
cold dust Aj 


4.0 idm" 

91.3 K, 1.92 Lo* 

4.2 mag * 


Varying, Best Fit, 


, Faint State 


composite T,, R^, L, 
Aj stellar extinction 
hot dust T, L 
Aj hot dust 


3514 K, 4.27 /?o, 2.49 Lq 
9.86 mag 
562 K, 1.10 Lo 
Omag'^ 


Varying, Best Fit, 


Bright State 



composite T,, R,, L. 3525 K, 4.66 Rq, 3.01 Lq 

A J stellar extinction 10.1 mag 

hot dust r,L 579 K, 1.63 Lo 

Ay hot dust 0.57 mag 



a Extinction wavelength dependence adopted from Becklin 
et al. (1978); Mathis (1990). A transition-wavelength initial 
guess of 3.5 jum was used in the fit, but was allowed to float freely 
in the fit to the faint state photometry. The fit for the bright state 
photometry fixed the transition-wavelength to that of the best fit 
in the faint state. 

b When we allow the cold dust temperature and luminosity to 
vary in our SED fit to the bright state photometry, the effect on 
the best fit parameters is marginal: <1 K, <1% luminosity, ~0.1 
mag A/ extinction difference. Thus, we fix these values in the 
bright state fit to reduce the degrees of freedom. The A./ value 
for the cold dust is larger than the best fit hot dust extinction, 
which seems counterintuitive. This may relate to non-blackbody 
thermal dust grain emission, or the cold and hot dust being asso- 
ciated with different components of the visual binary, c Negative 
extinction magnitudes were not allowed in the fit. 
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